Heavy fermions in an optical lattice 



Michael Foss-Feig, 1 ' 2 Michael Hermele, 1 Victor Gurarie, 1 and Ana Maria Rey 1,2 ' 3 



o 



00 
(N 



Department of Physics, University of Colorado, Boulder, Colorado 80309, USA 
2 JILA, Boulder, Colorado 80309, USA 
3 NIST, Boulder, Colorado 80309, 
(Dated: July 30, 2010) 



USA 



We employ a mean-field theory to study ground-state properties and transport of a two- 
dimensional gas of ultracold alklaline-earth metal atoms governed by the Kondo Lattice Hamiltonian 
plus a parabolic confining potential. In a homogenous system this mean-field theory is believed to 
give a qualitatively correct description of heavy fermion metals and Kondo insulators: it reproduces 
the Kondo-like scaling of the quasiparticle mass in the former, and the same scaling of the excitation 
gap in the latter. In order to understand ground-state properties in a trap we extend this mean- field 
theory via local-density approximation. We find that the Kondo insulator gap manifests as a shell 
structure in the trapped density profile. In addition, a strong signature of the large Fermi surface 
expected for heavy fermion systems survives the confinement, and could be probed in time-of- flight 
experiments. From a full self-consistent diagonalization of the mean- field theory we are able to 
study dynamics in the trap. We find that the mass enhancement of quasiparticle excitations in the 
heavy Fermi liquid phase manifests as slowing of the dipole oscillations that result from a sudden 
displacement of the trap center. 
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I. INTRODUCTION 

Ultracold atoms in optical lattices offer a clean and 
controllable setting for the experimental investigation of 
condensed-matter Hamiltonians. A number of remark- 
able experiments along these lines [TJ [2], including the 
observation of the Mott insulator to superfluid transition 
of the Bose-Hubbard model using bosonic alkali-metal 
atoms [3, and more recently the metal to insulator tran- 
sition of the Fermi-Hubbard model using fermionic alkali- 
metal atoms |H [5] , have thoroughly established the im- 
portance and feasibility of such optical lattice emulations. 
Recently it has been proposed that several unique prop- 
erties of fermionic alkaline-earth-metal atoms (AEMAs) 
make them particularly well suited for the simulation 
of a broad class of condensed matter Hamiltonians de- 
scribing the interplay between internal (spin) and exter- 
nal (orbital) electronic degrees of freedom [6HTU] . These 
two-electron atoms are more complicated than their one- 
electron alkali-metal counterparts, however in the last 
few years a great deal of progress has been made in cool- 
ing a variety of bosonic and fermionic isotopes to quan- 
tum degeneracy [TTHT7] . Motivated by these develop- 
ments, here we discuss cold-atom probes of the Kondo 
Lattice Model (KLM). 

The KLM is a canonical model for studying the in- 
teraction of conduction (mobile) electrons with localized 
(immobile) spin-1/2 scattering centers [18] . In the most 
common presentation of the model there is one localized 
spin on each lattice site, and the only interaction consid- 
ered is an on-site Heisenberg exchange between th e con - 
duction electrons and the localized spins [see Fig. l(a)| . 
In a recent paper [10 the implications of Ref. [6 for sim- 
ulating the KLM in one dimension (ID) were discussed. 
In particular, we suggested realistic experimental probes 



of physics occurring in several different regions of the 
ID phase diagram. Part of the analysis in [10] showed 
via a mean-field calculation that important properties of 
the KLM's paramagnetic phase are manifest in dipole os- 
cillations of the cold-atom system upon displacement of 
a confining potential. Here we present a detailed exten- 
sion of these mean- field calculations to a two-dimensional 
(2D) system, where the model is more closely related to 
real solid-state systems. 
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FIG. 1: (a) In the Kondo Lattice Model the conduction elec- 
trons (red) can hop from site to site, and they in tera ct with lo- 
calized spins (blue) via a Heisenberg exchange, (b) 2D mean- 



field ground-state phase diagram as constructed in [19], n g 
being conduction electron density and v being a dimensionless 
measure of the interaction strength. PM is a paramagnetic 
phase, in which the heavy Fermi liquid behavior is expected. 
The FM (AFM) phase is where RKKY interactions [20H22] 
generate ferromagnetic (antiferromagnetic) order among the 
localized spins. Exactly at n g = 1 the AF phase gives way to 
a non-magnetic insulating state for sufficiently large v [23 . 

For different values of the exchange coupling the KLM 
gives rise to very different physics. We focus exclusively 
on the case of antiferromagnetic exchange (favoring anti- 
alignment of the localized spin and a conduction electron 



on the same site), which has been studied for over three 
decades as a model for inter-metallic compounds with 
anomalously massive quasiparticle excitations, termed 
heavy fermion materials [18, 24, 25 . Despite the wealth 
of theoretical machinery developed for and brought to 
bear on the antiferromagnetic KLM, the structure of the 
ground-state phase diagram (in D > 1) is not yet well un- 
derstood. Neverthele ss, th e following qualitative picture 
has emerged [see Fig. |l(b)| . At sufficiently weak coupling 
the conduction electrons mediate a long-range magnetic 
interaction (RKKY) between the localized spins [2QH22] . 
The system will have magnetic order with details depend- 
ing on the lattice structure and the conduction electron 
density. At stronger coupling, the conduction electrons 
screen the localized spins by binding to them in sin- 
glets. The system is paramagnetic, and the excitations 
are thought to belong to a heavy Fermi liquid (HFL). 

A major outstanding challenge in KLM research is to 
understand the quantum phase transition that separates 
these two regimes, and this challenge could in principle 
be addressed by careful analysis of an optical lattice em- 
ulator. But before there can be any prospect for the 
exploration of unknown aspects of the phase diagram, 
it is imperative that we understand the manifestation of 
relatively well understood KLM physics in an inhomo- 
geneous setting — the ubiquitous cold- atom trap. With 
this goal in mind, we present an extension to the trap 
of a well known mean-field theory (MFT) appropriate 
for describing the HFL [19]. Some ground-state proper- 
ties can be determined with relatively little effort within 
local-density approximation (LDA). Already at this level 
important Kondo physics can be distinguished in the 
trap. For example, the Kondo insulator phase induces 
a shell structure in the trapped density distribution, and 
the emergence of heavy quasiparticles shifts conduction 
atom weight out to the so called "large Fermi surface". 
Outside of the range of validity of the LDA, for example 
in non-equilibrium situations, a numerical self-consistent 
diagonalization of the MFT is employed. Using the exact 
mean-field wave function we explore center-of-mass oscil- 
lations resulting from a sudden displacement of the trap 
center, and find a clear signature of the heavy fermion 
mass enhancement. 

The paper is organized as follows. In Section [IT] we 
briefly review technical aspects of simulating the KLM 
with AEMAs. In Section [m] we introduce the MFT 
and review its relevance to heavy fermion physics in a 
translationally invariant system. Section [TV| exploits the 
translationally invariant theory to characterize important 
ground-state properties of the trapped system within the 
LDA, including density and quasi- momentum distribu- 
tions. In Section |V| we explain how to implement self- 
consistent diagonalization of the inhomogeneous MFT, 
and in SectionlVllwe apply these solutions to calculations 



II. SIMULATION OF THE KLM WITH 
ALKALINE-EARTH-METAL ATOMS 

There are two primary features of AEMAs that make 
them suitable for simulation of the KLM: (1) The long- 
lived 1 6'o (g) and 3 Pq (e) clock states can be trapped in- 
dependently by two different but spatially commensurate 
optical lattices [26]. This enables us to consider localiz- 
ing one clock state in a deep lattice while keeping the 
other state mobile in a shallow lattice. (2) Both clock 
states have total electronic angular momentum J = 0. 
As a result there is essentially no coupling between the 
different nuclear spin states and therefore no spin chang- 
ing collisions (in contrast to alkali metal atoms). This 
enables us to consider an ensemble where just two states 
in the hyperfine manifold are populated, since to a very 
good approximation no transfer of states into the rest of 
the manifold can occur. We will identify these two pop- 
ulated levels with spin up and spin down of the electrons 
in the KLM. 

In order to make the above discussion more quantita- 
tive, we must introduce some notation of Ref. [6 to de- 
scribe and compare various energy scales. If the lattice 
containing the atomic state a G {#, e} has a potential V a 
and lowest Wannier orbital w ai then we define hopping 
energies J a = Jd D rw a (i-* - r) [-|^V 2 + V a ]w a (yj - r) 
(m is the atomic mass, and (r^, r^) are the centers of two 
nearest neighbor sites). In addition to these hopping en- 
ergies there are a variety of possible interactions, which 
at sufficiently low temperatures can be parameterized in 
terms of the 4 s-wave scattering lengths a ee , CL gg and af g , 
for scattering in the states |ee), \gg), and -4= (\eg) db \ge)) 
respectively. If we neglect interactions between atoms 
that are not on the same site, the four possible interac- 
tion energies are: U aa = (^h 2 a aa /rn) J d D rw^ (r) and 
U% = {^tfa%/m) f d D vwl (r) w* g (r). 

We are interested in the case where U ee ^> J e and 
U gg <C J g . Under these conditions, and with a suffi- 
ciently weak confining potential £1 (z), the e atoms form 
a unit filled Mott insulator at the center of the trap and 
the g atoms are to a good approximation non-interacting. 
We choose the e atoms to represent the localized spins be- 
cause they would otherwise suffer lossy collisions |B]. At 
temperatures well below U eei to lowest order we can drop 
both U ee and J e and simply constrain the Hilbert space 
to have one e atom per site. Certainly to higher order in 
J e /U ee there would be super-exchange interactions, but 
we assume these to be negligible compared to all other 
terms in the Hamiltonian. Defining V ex = (U+ — U~) /2 
and neglecting terms which are constant in consideration 
of the constraint on the e atom density, what remains is 
an inhomogeneous version of the Kondo Lattice Hamil- 
tonian: 

^K = ~Jg Yl C lg* C jg* + V exYl 
(i,j),cr %<j<j' 
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of non-equilibrium dynamics. Section VII concludes with 
experimental considerations. 
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In Eq (B n ia = £ 
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, where c] ao . creates an atom The resulting quadratic Hamiltonian is 



in the lowest Wannier orbital on the site centered at r^, 
in electronic state a E {e,g}, and (nuclear) spin state 
a. In the rest of the paper we will consider only D- 
dimensional hypercubic lattices, with a particular focus 
on D = 2. In the first term of the Hamiltonian we have 
used the convention that (i,j) restricts the summation 
to nearest-neighbor sites. It is important to observe that 
satisfaction of the inequalities U ee ^> J e and U gg <C J g 
places no fundamental constraint on V ex /J g . This inde- 
pendence of parameters is a unique feature of the AEMA 
simulations proposed in Ref. [6 . Previous proposals to 
study the KLM using alkali-metal atoms lack this tun- 
ability because they either populate multiple bands of a 
single lattice (Ref. [27], in which case U ee , U gg and V ex 
all scale with the lattice depth) , or generate Eq (IT]) as an 
effective Hamiltonian which necessarily operates at weak 
coupling (| V^x | / Jg <C 1) [28]. In the present scheme, V ex 
will vary from one isotope to another, and can be fur- 
ther adjusted by offsetting the two lattices with respect 
to each other (to decrease the overlap between Wannier 
orbitals). 



III. MEAN-FIELD THEORY 

There is no nontrivial parameter regime in which Eq. 
(fil) can be solved exactly, and approximations have to be 
made. In the HFL phase a qualitative understanding can 
be gained by treating the interaction term at mean-field 
level [19]: 

Vex c i ga c i eo -r c iga' c iea ~~ ** "ex / J Vj \ C i Q(J ' c ie(j' ■ C iea /C igcr / 
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where we have defined the on-site hybridizations 



Vi 
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This expectation value, along with all others in this pa- 
per, is taken in the ground state of the mean-field Hamil- 
tonian [see Eq. Q below]. Because the mean- field 
Hamiltonian is quadratic, the ground state is a Slater 
determinant of the lowest (N e + N g )/2 single particle 
states (half the total number of atoms due to spin de- 
generacy), where N a is the total number of a atoms. 
This decoupling is not unique, for instance there could 
be terms that mix different spin states, but if the KLM is 
generalized to allow the spin index a to run from 1 to N 
(N = 2 for the KLM), then the above decoupling is exact 
in the limit N — >• oo [29, 30 . One should keep in mind, 
however, that the decision not to keep any spin mixing 
terms in the decoupling guarantees that this MFT can- 
not capture a transition to a magnetically ordered state, 
and therefore cannot remain valid for arbitrarily small 
v (where RKKY-induced magnetic order should exist). 
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where the \±{ e are local chemical-potentials needed to en- 
force the constraint of one e atom per site on average 
(in the exact KLM we have the much stronger constraint 
((hie — I) 2 ) = 0> but in the MFT this must be relaxed to 
riie = (hi e } = 1). The theory is now paramagnetic, be- 
cause 7-Lm does not couple different spin states, and the 
eigenstates are created by quasiparticle operators: 



vt 
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v l c ■ 



(5) 



Here q E {1,2,...} is an index that labels the different 
quasiparticle eigenstates, and not a wave vector. Solving 
the MFT then means finding the above mode coefficients 
u l q and v l , but this needs to be done self-consistently: 

As parameters in the Hamiltonian the Vi determine the 
mode coefficients, but in turn the mode coefficients de- 
termine the Vi via the definition 
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(6) 



We sum over a number of modes equal to half of the 
total number of particles, but put two particles in each 
mode to account for the degeneracy of spin. Similarly, 
the onsite g and e atom densities are given by 
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q=l 

(N g +N e )/2 
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(7) 



If we define the ground-state energy E = (Km), the self- 
consistency condition and the local constraints r^ e = 1 
can be written compactly as 



dE/dVi = and dE/dfi ie = 



(8) 



respectively. 



A. Translationally Invariant Case 

Before solving the MFT in the presence of a trap, it 
is useful to review known results for the translationally 



invariant case (Q = 0) [19 . It is common to assume that 
the hybridizations and chemical potentials retain the dis- 
crete translational invariance of the exact KLM Hamil- 
tonian (Vi,fii e —> V,fi e ), in which case the Hamiltonian 
in Eq. (El) reduces to 



^ = -j g J2 4 



igcr c jgo- ' ^ e 



(ij)c 




+ V ex J2 V {4 g *Cie. + cLcigJ) - W ex V 2 , (9) 

J\f being the total number of lattice sites. It is common 
to include a g atom chemical potential by adding a term 
—fig ^2 i fiig in Eq. [9j but we find it conceptually simpler 
to work at a fixed n g for the time being (we will trade in 
rig for fig in the next section to apply the LDA). Due to 
the translational invariance the nonconstant part of Ht 
(those terms containing operators) can easily be diago- 
nalized for general V and fi e by going to Fourier space. 
Defining Ck aa = -h= V . Cj aa e lVj ' k (where Vj is the posi- 
tion of site j and k is a quasi-momentum vector in the 
first Brillouin zone), we find quasiparticles 



<*h* = C kga U ( k ) + 

and the quasiparticle spectrum 
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£±(k,V,fjb e ) = 



e(fc)+/i e 1 
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(11) 
Here e (ft) — — 2 J g Y^X=i cos(^a) (a being the lattice 
spacing) is the tight-binding dispersion for the g atoms 
in a D-dimensional hyper-cubic lattice, and the two 
branches of the quasiparticle spectrum are called the up- 
per and lower hybridized bands. If J\f is large then the 
energy per site of the translationally invariant system, 
given by E(V,fi e ,n g ) = (Ht) /A/", can be written as an 
integral over the Fermi volume T . For example, onaD 
dimensional lattice with n g < 1 the total energy per site 
is: 



E{V,fi e ,n g ) = -V ex V 2 -fi e 






k£-(k,V,fJL e ), 



(12) 

with the dependence on n g hidden in the limits of in- 
tegration. Only the lower hybridized band is populated 
because n g = 1 corresponds to completely filling the first 
Brillouin zone (since there is also one localized atom per 
site). For 1 < n g < 2 the lower hybridized band will 
be completely full and the upper band will be partially 
integrated over. For a given n g , the correct fi e and self 
consistent V can be found by solving the coupled equa- 
tions 



dE(V,fjb e ,n g ) 
dV 



and 



dE(V,fjb e ,n g ) 

8fl e 



0. 



(13) 



It will be useful for later application of the LDA to ap- 
preciate that the only knob we can turn to determine the 



MFT solutions (beyond the couplings J g and V ex in the 
Hamiltonian) is a choice of the g atom density. Once n g 
is chosen, the correct V and fi e are fixed by Eq (13). 
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Spectrum of the translationally invariant MFT for 
shown in ID for clarity. 



(a) and for 
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(b) 



The blue dashed line is the e atom band {E — fie), and the 
red solid curve is the g atom band (E — —2J g cos(fc)). The 
purple dotted lines are the hybridized bands S± , with the 
thicker solid section covering the Fermi volume. In |(b)| the 
lower band is filled and separated from the upper band by 
the hybridization gap Ah (green slashed region) , causing the 
system to be insulating for n g — 1. 



Defining the dimensionless coupling 

V = ~2V ex /Jg, 



(14) 



one finds that for small v the hybridization V will also 
be small. The two hybridized bands track the e and g 
atom bands closely, however mixing due to finite V causes 
the crossing to be avoided [Fig. ([2|]. For less than unit 
filling of the g atoms, the ground state is obtained by 
filling the lower hybridized band with quasiparticles out 
to the Fermi surface (kp = f (1 + n g ) in ID). The part 
of the lower hybridized band running nearly parallel to 
the line E — fi e [blue line in Fig. |2(a)| holds most of the 
e atom weight, which indicates that for decreasing n g we 
must pull fi e down below the center of the g atom band 
in order to satisfy n e — 1. This pins the Fermi surface 
to the flattened part of the hybridized band (called the 
Kondo resonance), which is the origin of the quasiparticle 
mass enhancement. 

When n g — 1 the lower hybridized band is completely 
full, and fi e = because of a particle hole sym metry [3T] . 
The system is an insulator with a gap Ah [Fig. |2(b)]. De- 



spite having the appearance of a band gap in the MFT, 



this gap should be understood as arising from correlation 
effects: Filling the lower hybridized band means having 
one conduction atom per unit cell, so non-interacting 
band theory of the conduction atoms would predict a 
metallic state. From Eqs. Q and (|6| one can see that 
at v = oo and fi e = we have V = 1. It then fol- 
lows from Eq (11) that the MFT predicts Ah ~ J g v at 



strong couplingTlYt weak coupling, more careful analysis 
shows that the MFT predicts a non analytic gap scaling 
Ah ~ J g e~ x l v . This weak coupling expression should be 
understood as correct for v small but not less than some 
critical v c > 0. The critical coupling v c marks the phase 
transition between the HFL and a (presumably) magnet- 
ically ordered phase, in which the MFT is not valid [see 

Fig. [T(b)1. 

It has been implicit in the discussion above that both 
the e and g atoms are to be included in the Fermi volume 
of the MFT ground state. However, one should keep in 
mind that the actual size of the Fermi surface is a subtle 
and important issue in the KLM [3TJ [32] . Looking at Eq. 
([I]) , one might expect the volume of the Fermi sea in any 
KLM ground state to be determined by the number of g 
atoms, since the e atoms represent localized spins with no 
charge degrees of freedom. However, the KLM at small v 
is an effective model for the Periodic Anderson Model at 
large Z7, for which Luttinger's theorem [33] is expected 
to hold [32] (meaning that the Fermi volume should be 
determined by the number of g and e atoms). These two 
scenarios are referred to respectively as a small and large 
Fermi surface. The HFL ground state, which is described 
by our MFT, is known to have a large Fermi surface. The 
Kondo insulator is a cousin of the HFL, which arises when 
the density of g atoms increases to unity, and the large 
Fermi surface of the HFL completely fills the Brillouin 
zone. Our focus in this paper is on the HFL and Kondo 
insulator, largely because we believe it is likely that the 
HFL and Kondo insulator phases will occupy a significant 
part of the param agnetic region of the phase diagram 
[PM in Fig. |l(b)| . Indeed, on the square lattice that 
we consider the Kondo insulator is known to occur for 
v > 1.45 (when n g = 1) [23]. Given the close relationship 
between HFL and Kondo insulator, the HFL most likely 
at least exists nearby (i.e. for v > 1.45 and n g < 1). 



IV. LOCAL-DENSITY APPROXIMATION 

The simplest way to extend the MFT to an inhomo- 
geneous system is via the local density approximation 
(LDA). The idea is that on a given site, we consider the 
trapping potential to be a chemical potential in a trans- 
lationally invariant system with Hamiltonian 



^T-HgJ2 c l 
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(15) 



with /jijg = fi g — Q(j). On physical grounds one might 
argue that the local chemical potential should be applied 
to the e atoms as well. However this is unnecessary; the 



e atom density is fixed to be one per site, so including 
the effect of the trap on them simply provides an overall 
constant energy shift (more technically, if it were applied 
to the localized atoms as well the local jij e would simply 
readjust to absorb it). We now have one translation- 
ally invariant problem per lattice site, each of which can 
be solved as described in Section IIIII There is a minor 



complication owing to the choice in Section III to work 
at fixed density, not fixed chemical potential. But this 
can be resolved easily by remembering that at zero tem- 
perature the chemical potential is the energy required to 
add one quasiparticle at the Fermi surface. Solving the 
translationally invariant model at fixed fij g simply means 
finding the rij g satisfying 



dE(V,fi e ,n jg ) 



39) _ 



dUjg 



^jgi 



(16) 



since the left hand side is the energy cost of adding one 
quasiparticle at the Fermi surface. The derivative on the 
left hand side can be expanded as 



dE dV l dE dn e 
dV dn« 



dE 



dfi e drijg drij 



(17) 



and the first two terms are zero for the MFT solution by 
Eq. (13). The last term is the change in the energy, at 



fixed V and fi e , due to the addition of a single quasipar- 
ticle at the Fermi surface, so we conclude that 



fi jg = £±(k F (n jg ),V(n j g),fi e (n jg )) 



(18) 



Here kp is any vector falling along the Fermi surface con- 
sistent with the filling n^, and the ± means to choose 
the branch that yields a solution: there will only be one, 
since both branches are monotonic functions of n g and 

H 39 



they are separated by a gap. Solving Eq. (18) for n 



on each lattice site, while adjusting fi g to obtain the cor- 
rect total number of particles, constitutes the LDA so- 
lution. The LDA gives us immediate access to many 
important groundstate properties, such as the real-space 
and momentum-space density distributions. 



A. Real space g atom density distribution 

In order to understand the qualitative features of 
trapped g atom density distributions it is helpful to plot 
the ground-state mean-field phase diagram in the fi g —v 
plane [Fig. [3J. Moving from the outside of the g atom 
cloud towards the center of the trap corresponds to mov- 
ing from the bottom to the top of the phase diagram 
along a line of fixed v. For sufficiently small N g , we will 
reach the trap center before entering the region of Kondo 
insulator, and therefore we expect to be everywhere in 
the rig < 1 heavy fermion metallic state. If, however, we 
choose N g sufficiently large, we will breach the n g — 1 
Kondo insulator phase and a unit filling plateau should 
develop at the center of the trap. Increasing N g further 
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V 
FIG. 3: Schematic ground- state mean- field phase diagram 
of the KLM as a function of chemical potential fi g and di- 
mensionless coupling v. The black dotted lines are given by 
li g = ±Ah/2, which crosses over from Kondo like scaling 
(Ah ~ J g e~ x l v ) at small v to linear scaling (Ah = J g v) at 
2;e v. 



we will find that near the center of the trap we exit the 
Kondo insulator phase into the n g > 1 heavy fermion 
metallic phase. This behavior is very similar to the Mott 
plateaus of the Hubbard model [34j [35], only here we 
have a Kondo insulator and not a Mott insulator. We 
also emphasize an unusual feature of the shell structure, 
that it exists in the density distribution of g atoms, which 
do not interact with each other directly. 
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FIG. 4: (a) The formation of a Kondo plateau be twee n v\ 
and r<i at high trap fillings, obtained within LDA. (b) The 
hybridization also displays the plateau, but as the density 
incre ases pas t n g — 1 the hybridization goes back down. In 
both |(a)| and |(b)| the parameters used are N g w 550, Q/J g — 
45/1000, and v = 8. 



In Fig. |4(a)| we demonstrate such a plateau for a sym- 
metric trap geometry Q(i) = fir 2 (where r is the distance 
of site i from the origin, in lattice units). This result was 
obtained by exploiting the rotational symmetry and solv- 
ing a simpler ID problem in the radial coordinate r. We 
are still solving Eq ( |l8| ), but now instead of fijg we use 
li rg = fig — fir 2 , andwe rotate the solution around the 
axis of cylindrical symmetry to obtain the plot in Fig. 
4(a) Of course for this to work the g atom cloud radius 



must be relatively large in lattice units, otherwise the 
lack of rotational symmetry of the lattice itself would be 
relevant. Th e plat eau also shows up in the hybridization 
profile [Fig. |4(b)| , however near the center of the trap 
the hybridization dips back down. Physically this occurs 
because for n g > 1 some g atoms are forced to pair with 
each other into singlets, reducing the number available 



to hybridize with the localized atoms. 

Defining 7*1(2) to be the interior (exterior) radius of 
the plateau, it is clear from the argument used to link 
Fig. [3] to Fig. [4] that ft (r| - r 2 ) = A H . More gener- 
ally, in an exact treatment of the KLM the LDA solution 
will predict ft [r\ — r 2 ) = 2A qp , where A qp , called the 
quasiparticle gap, is the difference in energy between the 
unit-filled ground state and the lowest energy state with 
one g atom added or removed (these are both the same 
because of particle hole symmetry) [31] . Physically this 
result reflects the LDA assumption that transfer of den- 
sity from one site to another does not affect the energy to 
first order (dE(nj g )/drij g + fl(j) = H g is constant across 
the trap). Therefore the energy needed to add a quasi- 
particle at 7*1 (A qp — n g + fir 2 ) must be opposite to the 
energy needed to remove one at r2 (A qp + fig — fir 2 ), 
giving ft (r| - r 2 ) = 2A qp . 

At strong coupling, where the plateau is most visible, 
the exact quasiparticle gap of the KLM can be calculated 
by considering a single site (due to the relatively small 
hopping) [31 . If on one site we have n g = n e = 1 initially, 
adding or subtracting a g atom costs (3/4) J g v = A qp . 
Therefore we can be certain that for large v the plateau 
size should satisfy ft(r 2 — r 2 ) = (3/2) J g v. In this limit 
Ah = J g v, so the MFT underestimates the size of the 
plateau in the strong-coupling limit. 



B. g atom quasi- momentum distributions 

At mea-field level a large Fermi surface arrises by as- 
sumption: By assigning a nonzero value to the hybridiza- 
tion matrix element V, the e atoms are liberated into the 
Fermi sea. This is well known, and proves nothing about 
the Fermi surface in the actual KLM ground state. How- 
ever, as an example of how the structure of the large 
Fermi surface survives in the trap, we proceed to calcu- 
late the quasi- momentum distribution in the LDA [36J. 
In the translationally invariant MFT, there is a unique g 
atom quasi-momentum distribution 



h[k,n g ] = f[k,n g ] \u{k)\ z 



(19) 



associated with every possible filling n gi where /[fe, n g ] is 
the zero temperature Fermi function for non-interacting 
electrons at a filling 1 + n g . The LDA approximation 
to the quasi-momentum distribution is then obtained by 
summing n over the densities on each lattice site 



n L DA[k] = y^n[k,n jg ], 



(20) 



where the Uj g are calculated as discussed in the previous 
section. 

The existence of the trap introduces an unavoidable 
ambiguity in the definition of the small and large Fermi 
surfaces, since the filling is affected by what one chooses 
as the system size (i.e. should the sites outside of the 
g atom cloud be included?). However we find it reason- 
able, for the sake of comparison, to choose the system 
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FIG. 5: (a) The first Brillouin zone of a 2D square lattice. The central white region (within the dotted circle) is the small Fermi 
sea of the translationally invariant KLM with a conduction atom density n g ~ 0.19. The shaded green region (within the solid 
red line) fills out the large Fermi sea (the red solid line is the large Fermi surface, as defined in Section IV B). The red dotted 
line (just outside of the red solid line) is the large Fermi surface corresponding to a uniform filling equal to the density at the 
center of the trap. Even when there are no conduction atoms the large Fermi sea occupies half of the zone (within the dotted 
diamond), since there is still one localized atom per site contributing to the volume, (b) Quasi-momentum distribution of the 
conduction atoms plotted over one quadrant of the Brillouin zone, showing a prominent feature at the large Fermi surface (red 
ribbon), (c) A ID cut along the line k y = k x in (b), hLDA[k x , k x ] is the blue solid line, and the large Fermi surface is the red 
dotted line. For all plots the parameters used are N g w 116, Q/J g = 5/1000, and v = 8. 



size to be defined by the diameter of the g atom cloud. 
Then the trap can be considered as a perturbation on the 
translationally invariant system that is sufficiently small 
as to not force the density to zero anywhere. The large 
Fermi surface then belongs to a translationally invariant 
system of this size containing the same number of con- 
duction atoms as there are in the trap. Applying this 
definition to ~ 116 conduction atoms in the symmetric 
trap of Section 



IV A 



(this time with Q = 5J^/1000), we 
plot the large Fermi surface as a solid red line in Fig. 
[5|a). In Fig. [5jb) we plot the quasi-momentum distribu- 
tion riLDA[k] for the trapped system alongside the large 
Fermi surface of the translationally invariant system. 

It is worth appreciating that the large Fermi surface ac- 
tually survives smearing by the trap better than it would 
for free fermions. This is in part because in the averaging 
over the trap only Fermi surfaces between the black and 
red dotted lines of Fig. [5Va) contribute (for free fermions 
the averaging would contain arbitrarily small Fermi sur- 
faces). Another reason is that, as we approach the edge 
of the conduction atom cloud and the local Fermi surface 



7r, the height of the 



recedes toward the line \k x \ + \k 
discontinuity also decreases to zero. 



V. SELF-CONSISTENT DIAGONALIZATION 

In order to probe mass enhancement in the heavy 
fermion state we would like to study dynamics of the 
conduction atom cloud in response to a trap displace- 
ment. This type of dynamics has been shown to be an 
excellent diagnostic tool for strongly correlated systems 
in experiments using bosonic as well as fermionic atoms 
At mean-field level this requires us to know 



the self-consistent ground state; we must move beyond 
the LDA and self-consistently diagonalize Hm- Analytic 
progress is made difficult by the necessity to maintain site 
dependent hybridizations, and we proceed numerically. 
We start by making an initial guess for the Vi and /i^ e 
using LDA. We then diagonalize Hm numerically, calcu- 
late new Vi from Eq. (|3| , and repeat until the Vi converge 
to a self consistent value. This process is complicated by 
the need to repeatedly determine the /i^ e that satisfy the 
constraint of one e atom per site on average (these will 
change every time the Vi are updated). A discussion of 
the procedure we use to determine the correct \ii e can be 
found in the Appendix. 

To simplify the numerics, we consider a geometry 
where the potential changes only in the longitudinal di- 
rection (length aLi and g atom hopping J/), while the 
transverse direction (length aL t and g atom hopping 
J t ) is homogenous with periodic boundary conditions. 
The motivation for this unusual geometry is to retain 
a number of constraint equations equal to the longitudi- 
nal length, thereby accessing relatively large system sizes 
(up to ~ 1000 sites is reasonable). Adopting a notation 
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ci„„„„ for the creation operators, where m labels the site 



in the longitudinal direction and n labels the site in the 
transverse direction, we define partially Fourier decom- 
posed operators 



J 



1 

7r t 



£< 



(21) 



Rewriting the nonconstant part of the mea- field Hamil- 
tonian [Eq. Q] in terms of these new operators, it 
decouples into the sum of L t effectively ID Hamiltoni- 
ans, which we label by the transverse quasi-momentum 



y L t a^ ' a j 



£"* 






-^ E4««,w+^«EM C * 



kga ikea 



A r — V 

^ikea^ikga v i 



/i) + J2 (Mi 



1) ct c 



ikea ikea 



+ £(-2J t cos(fca) + fit 2 )4 flff c ifcfl<T . (22) 

One can give a simple physical interpretation to these 
new Hamiltonians: the energy gained by a g atom dis- 
persing with quasi-momentum k in the transverse direc- 
tion [2Jt cos(fca)] has been incorporated as a chemical 
potential in a ID model describing hopping in the longi- 
tudinal direction. The eigenstates are created by quasi- 
particle operators 



qka 






y U qk C ikga 



' J qk 



ikea 



(23) 



Similar to Eq. ([5|, here q E {1, 2, ...} is an index labeling 
the eigenstates in order of increasing energy (for a given 
k). Unless the transverse bandwidth is smaller than the 
level spacing of the harmonic potential (4J t < 2y/QJi), 
the different k modes cannot be expected to have the 
same occupation numbers, which leads us to define k- 
dependent Fermi levels q^. For a given k, q^ is the small- 
est q for which the mode labeled by q and k is unoccupied. 
The ground state is 

i*)=nn4i°>. ( 24 ) 

ko- q<q% 

and the expectation value in Eq. ([3| becomes 

^ = EEka+w ( 25 ) 

k q<q k F 

VI. DYNAMICS 
A. Single particle dynamics in a trap 

Before introducing dynamics in the MFT, we briefly 
mention some important details of single-particle dynam- 
ics in a lattice plus harmonic potential [40-43 . If there 
are no interactions, the trapping potential is parabolic 
with curvature Q, and we work in ID for simplicity, then 
the g atoms are described by the Hamiltonian 



Ho 



-Jo 



E 



^E 



1 C iga C iga' 



(26) 



riodic Mathieu functions, with the corresponding eigen- 
value E n determined by the function's characteristic pa- 
rameter [40]. The solutions are complicated, but in cer- 
tain limits they have a very simple form. We will ex- 
clusively consider the limit where q = 4J g /Q ^> 1. Then 
1/y/q is a natural small parameter for expanding both the 
eigenfunctions and eigenvalues, and to lowest non-trivial 
order one obtains: 



c 



E r ' 



V2 



2™n! (<?7T 2 ) 1/4 
■l-2q + 4V9 



exp 



(2n+l) 2 



8 



W) 



where £ = i (4/g) / and H n are Hermite polynomials. In 
this approximation the expectation value for the center 
of mass (COM) of the ground state, evolving due to a 
small displacement 8 (measured in lattice units), can be 
obtained in closed form (the motion of excited states is 
similar but slightly complicated by various combinatorial 
factors): 



(X(t)) = 5 exp 



s 2 fnty 



COS 



Ah) 1 



2x 2 



sin 



.(28) 



W4 : 



In the above x = (#/4) ^ is a characteristic oscillator 
length (in lattice units) and uj* = fl^/q/h is the charac- 
teristic frequency. It is interesting to note that there is 
periodic (non-dissipative) damping of the oscillation am- 
plitude. In the absence of the lattice this damping would 
be forbidden by the generalized Kohn theorem [44], ac- 
cording to which COM oscillations in a harmonic poten- 
tial should be undamped under fairly general conditions. 
It is important for our purpose in what follows 
to observe that the tunneling matrix element is in- 
versely proportional to the lattice-effective mass m* = 
(2J g a 2 /h 2 )~ 1 . Therefore, the effective oscillator period 
(at which, for sufficiently large q, all low lying sin- 
gle particle modes oscillate) is related to the mass as 
r * _ 2jl ^ yVn*. In what follows, we take the periodic- 
ity (f ) of oscillations to define the effective mass for the 
interacting system (m): m/m* = (f/r*) . 



B. MFT dynamics 

Once we have solved for the MFT groundstate, cal- 
culating dynamics is relatively straightforward. One can 
easily write down a time dependent Schrodinger equation 
for the mode coefficients. Just as ip q (x) — » ip q (x,t) in 



A single particle eigenstate |^) = ^^fcj |0) has continuum quantum mechanics, here we have a discrete 
coefficients ^f given by the Fourier components of pe- analogy u J qk — » u J qk (£), and likewise for the v J qk . The time 



dependent discrete Schrodinger Eq. reads 



ih- 



du{ 



qk 



dt 



Ji^ + u^-nu-sfu 



J 
qk 



2J t [cos(fea)] u 3 qk - V ex Vji 



qk 



-ih- 



dvi 



qk 



dt 



"fi (J ~ $f <k ~ Hev j qk - V ex V jU H29) 



Eq (29) can be obtained formally by considering the co- 



efficients in Eq (23) to carry the time dependence and 



then solving the Heisenberg equations of motion for the 
quasiparticles 



a— -a 



dt i k ° 



t _ 



Wfc.ajfc, 



(30) 



(H 5 k being a shift of %k by S lattice sites). 

There is one subtlety in the time evolution: although 
the exact KLM Hamiltonian has a local U(l) symmetry 
to conserve the e atom density distribution, the MFT 
does not: the occupation of the localized band will cer- 
tainly evolve in time, which we cannot allow. The solu- 
tion is to let the chemical potentials be time dependent. 
It is easy to check that the first time derivative of the e 
atom density does not depend on the chemical potentials, 
and hence the constraints drii e /dt = cannot determine 
the \ii e . Instead they must be chosen to satisfy the sec- 
ond order constraints d 2 ni e /dt 2 = 0. It follows directly 
from Eqs. (29) that dn^jdt = initially, since all the 



mode coefficients can be chosen to be real. Keeping the 
second derivative zero at all times means the first deriva- 
tive does not change; it will always be zero and the local 
constraints will be obeyed. To find an explicit formula for 
the chemical potentials we first express the constraints in 
terms of the mode coefficients: 



dt 2 '' 



= 2 J2 J2 {%k% k + v\ k % k + 2vj fc v* fc ) . (31) 



q<q% 



Using the Schrodinger equation to evaluate the time 
derivatives in Eq. (pll), some algebra leads to: 



Mi, 



-2v ex n [dif ( Ci - fi) - m [di K + h)] + n [d iXi 



K[di 



■[di. 



(32) 

where 5ft and ^ are the real and imaginary parts respec- 
tively. The latin letters are defined as: 



(Li = 


= Tl ^2 ^\k U qk 


bi = -kT.'*U u % 1 


di ~- 


= TT t zZ ^qk U qk 


1 / 

fi = T~ t ^ %k V qk 


Ci ~- 


= TT t S U l qk U qk 


*i = ~^E v qk ul k cos(ka) 



where £ =E k E q<q %- 



Using Eq. (|25[) and (32) to write the V\ and \±i e in 



terms of the mode coefficients reduces Eq ( 29 ) to a sys- 



tem of 4:L t Lf coupled first order differential equations, 
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FIG. 6: (a) Here we plot the time evolution of the COM after 
displacement of the trap by one lattice site (red line). The 
blue dashed line is the noninteracting dynamics for compar- 
ison. The parameters used in t his p lot are q — 235, v — 4, 
L t = 8, N g = 40, and Ji/Jt = 2. |(b)| The effective mass scales 

by extracti ng t he 
at 



like rh/m* ~ (j/r*) (see the vertical axis) 

period of oscillation from several traces like the one in (a) 



several different values of v we find significant enhancement 
of the effective mass as v ~-> 1 . The red circles are from self- 
consistent diagonalization of the MFT, and the blue line is a 
guide to the eye. 



which can be integrated with standard methods. We cal- 
culate the dynamics resulting from displacement of the 
trap by one lattice site for fixed J^Jt, and Vt at several 
different values of V^x^and extract the period of COM 
oscillations (see Fig. p|. Since rh/m* = (f/r*) , Fig. 
* n demonstrates the emergence of heavy quasiparticles 



|6(b)1 ( 



as we decrease v. 

In the thermodynamic limit of the translationally in- 
variant MFT, the density-of-states effective mass truly 
diverges as v —> 0. This divergence is not physical, be- 
cause for sufficiently small v RKKY interactions lead to 
the development of magnetic order, at which point the 
paramagnetic MFT is no longer valid [18 . In the trap, 
however, the effective mass does not diverge in the small 
v limit, even within the MFT. This follows from the 
discreteness of the spectrum in the non-interacting the- 
ory; when the hybridization energy VV ex becomes smaller 
than the energy spacing, it cannot significantly change 
the system properties. Instead, as we decrease v the rel- 
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FIG. 7: COM oscillations for decreasing values of v, show- 
ing how the fast free-particle dynamics emerge on top of the 
slow quasiparticle dynamics. The red solid line is from MFT 
dynamics, and the blue dotted line is the noninteracting so- 
lution. The parameters used in this plot are q = 235, L t = 1, 
N g = 6, and Ji/J t = 2. 



atively fast free particle oscillations emerge on top of the 
slow quasiparticle oscillations. Eventually, the dynamics 
converges to that of the noninteracting system, as shown 
in Fig. 



VII. EXPERIMENTAL CONSIDERATIONS 
AND OUTLOOK 



adiabatically tightened until n g > 1 at the trap center. 

The observation of dipole oscillations over a small num- 
ber of lattice sites (< 8) has precedent in several exper- 
iments done with alkali- metal atoms [37, 38 . Usually 
the COM velocity of the cloud as a function of time is 
measured directly by TOF imaging, and the COM posi- 
tion is inferred from simple kinematics. There are also 
promising proposals to measure the COM expectation 
value in a non-destructive and precise manner by cou- 
pling the atomic motion to unpumped cavity modes [45] 
(the COM motion is inferred indirectly by its relationship 
to measurable quadratures of the photon field.) 

For the above mentioned experiments to be feasible, it 
is important that the requisite temperatures are within 
the reach of current technology. In real metals heavy 
fermion behavior develops only below the Kondo tem- 
perature Tk ~ J g e~ x l v . For u«l, which is always the 
case in heavy fermion metals, the exponential suppres- 
sion guarantees that Tk is always well below the Fermi 
temperature of the non-interacting conduction electrons, 
and that the mass enhancement (which scales as T^ 1 ) 
is extremely large. However one should not think of the 
smallness of Tk as being fundamental to heavy fermion 
behavior, but simply as the result of the bare parameters 
(v <C 1) which happen to exist in real metals. A larger 
Kondo temperature can always be obtained by choosing 
a larger v, at the cost of reducing the mass enhancement. 
By artificially creating the KLM in an optical lattice and 
choosing v ~ 1 the Kondo temperature could be esti- 
mated as Tk ~ J g ~ V ex . Therefore the temperature 
need only be smaller than the interaction energy, which 
is a much better situation than one encounters in propos- 
als to observe super-exchange or RKKY type physics. 



In order to obtain the heavy fermion metallic ground- 
state, it is necessary to prepare a Mott insulator of e 
atoms coexisting with a dilute cloud of g atoms. To real- 
ize the SU(2) symmetry of the KLM under consideration 
only two hyperfine states should be initially populated. 
In consideration of the lossy nature of e-e scattering, the 
starting point we envision is a Mott insulator of g atoms 
with doubly occupied sites in the center of the trap. Such 
a configuration can be obtained by making the lattice for 
g atoms sufficiently deep and the parabolic trap suffi- 
ciently tight. The deep lattice for the e atoms is un- 
occupied but already established. Taking advantage of 
the mean-field energy shift associated with double occu- 
pied sites, it is in principle possible to make the transfer 
\gg) —} (| eg) + \ge)) /y/2 with high efficiency. The sites 
occupied by a single g atom (at the edge of the trap) can 
also be addressed independently, transferring the atoms 
on them into the e lattice. If the g lattice is adiabatically 
lowered, we are left with the configuration desired. For 
sufficiently large v this procedure can also yield a unit 
filling plateau at the trap center (the Kondo insulator), 
but to observe the secon d layer in the shell structure dis- 
cussed in Section IV A (which has 3 atoms per site in 
the center of the trap, 1 e and 2 g), the trap should be 



VIII. SUMMARY 

This paper is motivated by recent progress in cooling 
and controlling AEMAs [TTHT7] . and proposals to use 
cold AEMAs to implement optical lattice simulations of a 
broad class of Hamiltonians [6H3 HH] 5 of which the KLM 
is one example. Building on these prospects, we make 
clear several cold-atom experimental signatures which 
would unambiguously demonstrate the physics found in 
heavy fermion metals and the closely related Kondo in- 
sulators. 

Local-density approximation has been used to extend 
a previously studied mean-field theory to an inhomoge- 
neous setting. At this level of approximation it is clear 
that the g atom density distribution holds a smoking 
gun of the insulating behavior of the KLM at half fill- 
ing (in the form of a density plateau), and it is plausible 
that the emergence of a large Fermi surface will be ev- 
ident despite the inhomogeneity of the trap. A more 
careful self-consistent diagonalization of the MFT with 
a trapping potential in two dimensions is explained in 
detail, and yields the ground state MFT wave function. 
We then show how this ground state can be time evolved 
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self consistently upon displacement of the trapping po- 
tential, and find that heavy fermion mass enhancement 
manifests as a slowing of the dipole oscillations when the 
coupling v is decreased. 
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Appendix: satisfying the local constraints 

For the initialization of the self-consistent diagonaliza- 
tion procedure, and anytime during the procedure im- 
mediately after the Vi have been updated, it is necessary 
to determine the /i^ e that satisfy the local constraint on 
the e atom density (n e i = 1). Finding the correct \±{ e 
amounts to finding a root of the vector- valued function: 



(^ e ) = i-2j2J2 



U qk 



(A.1) 



q<q% 



since Aj is the deviation from 1 of the number of localized 
atoms at site j. The dependence of A on the chemical 
potentials is implicit in the v J k , because they are ob- 
tained by diagonalizing the Hk, in which the \ii e appear 
explicitly. We proceed to solve Aj (/ii e ) = via New- 
ton's method, where the gradient is calculated in first 



order perturbation theory (this is just linear response 
theory) [7] : 



dfl ie 



2 
Lt 



EE E^ 

q'>q% 



q<q% 



^4'^'^ 



, (A.2) 



!R being the real part. Newton's method amounts to 
making successive linear approximations to the nonlin- 
ear equation Aj (fii e ) = 0: we choose a starting point /i^ e 
from the LDA, and then solve 



a,M.)-$>,.(!^ 



(A.3) 



for the 5 [lie by inverting the matrix d/S.j/djii e . We then 
change ji Q ie — » /i? e + <5/^ e , and repeat until the constraints 
are satisfied to within a desired tolerance. 

This procedure works wherever the hybridizations Vi 
are finite, but breaks down, for instance, near the edge 
of the trap. In general, we initiate the self-consistent 
diagonalization with a trap strength that is sufficiently 
weak to allow finite g atom density at the trap edges, 
and then turn up the trap slowly throughout the iter- 
ation. Whenever the hybridization at the edges begins 
to vanish, the Fermi level is flanked by two nearly de- 
generate combinations of a single e atom at the leftmost 
or ri ghtmost lattice site. The gradient defined in Eq. 
( |A.2[ ) becomes badly behaved at this point, because the 
degenerate states split the Fermi level, and they must 
be excluded from the sum (one can avoid the resulting 
divergences by exploiting the reflection symmetry about 
the center of the trap, but inclusion of these degenerate 
states nevertheless increases the precision necessary to 
reliably invert the gradient matrix). 
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